Prerequisite: it is assumed in this tutorial that time-lapses have been processed using the TissueMiner automated workflow
We provide one demo dataset (demo) and three processed time-lapses corresponding to the distal part of 3 wild-type Drosophila pupal wings (WT_1, WT_2 and WT_3).
No programming skill required: the tutorial Learning_the_R_basics_for_TissueMiner will guide the user, without any previous knowledge in programming, through the main aspects of data manipulation and visualization using the R language.
The power of an existing grammar to manipulate and to visualize data. The open-source R language now provides very efficient and simple ways to handle important amount of tabular data, using a so called “grammar” of data manipulation (dplyr package) and visualization (ggplot2 package) (by Hadley Wickham). We found R to be ideal for querying our relational databases (RSQLite package), and to perform subsequent calculations and analyses for understanding tissue morphogenesis. Thus, the user keeps control of his data by learning a few verbs and the simple syntax of the grammar. This gives the advantage of full autonomy on the data that can manipulated without restriction.
TissueMiner extends this grammar by providing an Application Programming Interface (API) to facilitate the visualization and quantification of cell dynamics during tissue morphogenesis. This API consists of a series of easy-to-use commands that we describe here on examples.
TissueMiner is compatible with the Python programming language. We provide a Python-tutorial of TissueMiner.
High performance computing platform. TissueMiner has been designed to be used in command-line in order to easily batch repetitive tasks and to run tasks on a high performance computing platform (cluster).
A comprehensive HOWTO to visualize and quantify cell dynamics in 2D tissues, by a simple copy-paste of the code encapsulated in shaded grey boxes into a R-shell such as in RStudio. We also provide a single script analyze_movie.R to perform tasks of the section 2 of this User Manual. We strongly encourage the user to use RStudio or a similar software to conveniently run the R programs step by step.
TissueMiner can easily be extended by the user to address project-specific questions
Tutorial layout: the code is encapsulated in shaded grey boxes that delimit code chunks. The results are displayed immediately below in an open box in which all lines start with the ## sign. Within code chunks, comments are visible in green and are prefixed with at least a # sign. The comments indicate what the following code does. Here is an example:
# This is a comment: the code below will print "Welcome to TissueMiner"
print("Welcome to TissueMiner")
## [1] "Welcome to TissueMiner"
sm allsm -nsm make_dbsm shear_calculate topo_countt1 polygon_class tri_categorizetm alias for simplicity. Example:
alias tm='docker run --rm -ti -v $(dirname $PWD):/movies -w /movies/$(basename $PWD) etournay/tissue_miner'tm sm shear_calculate topo_countt1 polygon_class tri_categorizesm command) builds the database and generates additional data about cell lineages, cell topology, cell geometry, and cell contributions to tissue deformation./usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
brew update
brew install git
export TM_HOME="${HOME}/tissue_miner"
git clone https://github.com/mpicbg-scicomp/tissue_miner.git ${TM_HOME}
${TM_HOME}/installation/osx/install_tm_api.sh
docker run -d -p 8787:8787 -e USER=rstudio -e PASSWORD=rstudio -v $HOME:/home/rstudio/share etournay/tissue_miner
-v option. Example:-v $HOME/my_R_scripts:/home/rstudio/scripts -v /media/fileserver/movie_repository:/home/rstudio/movies
echo "http://$(docker-machine ip default):8787"
algnModel variable# To define the location of TissueMiner scripts and configuration file, just open a terminal and type in:
echo "TM_HOME=~/tissue_miner" >> ~/.Renviron
echo "TM_CONFIG=~/tissue_miner/config/flywing_tm_config.R" >> ~/.Renviron
# On Windows, TM_HOME and TM_CONFIG must be defined in the .Renviron file.
# Jut type in: help("Startup") in a R session (Rstudio), and follow the manual
# To define the location of TissueMiner scripts and configuration file, just open a terminal and type in:
echo "TM_HOME=~/tissue_miner" >> ~/.Renviron
echo "TM_CONFIG=~/tissue_miner/config/default_config.R" >> ~/.Renviron
# To define the location of TissueMiner scripts and configuration file, just open a terminal and type in:
echo "TM_HOME=~/tissue_miner" >> ~/.Renviron
# Here provide the path to your configuration file
echo "TM_CONFIG=path_to_my_config_file" >> ~/.Renviron
# Define path to all processed movies: MUST BE EDITED BY THE USER
movieDbBaseDir="/Users/retourna/example_data"
# Define path a particular time-lapse called "demo"
movieDir <- file.path(movieDbBaseDir, c("demo"))
# Define a working directory where to save the analysis
outDataBaseDir=file.path(movieDir, "output_analysis")
# Set up path to the TissueMiner code
# This command requires that the global environment TM_HOME is defined in the .bash_profile
scriptsDir=Sys.getenv("TM_HOME")
# Load TissueMiner libraries
source(file.path(scriptsDir, "commons/TMCommons.R"))
source(file.path(scriptsDir, "commons/BaseQueryFunctions.R"))
source(file.path(scriptsDir, "commons/TimeFunctions.R"))
# Overwrite any default configuration for running the examples of the User Manual
source(file.path(scriptsDir, "config/flywing_tm_config.R"))
# Set up working directory
mcdir(outDataBaseDir)
| Functions | Description | Project |
|---|---|---|
| print_head | head the current table with the row number | TissueMiner |
| dt.merge | fast merging of two dataframes, possibility to suffix column names | TissueMiner |
| openMovieDb | open a connection to a movie database | TissueMiner |
| multi_db_query | aggregate data into one dataframe | TissueMiner |
| coarseGrid | assign grid elements to cell positions | TissueMiner |
| smooth_tissue | average quantities in time using a moving window | TissueMiner |
| align_movie_start | align movies at earliest common developmental time | TissueMiner |
| chunk_time_into_intervals | undersample time for local time averaging | TissueMiner |
| synchronize_frames | find closest frame to user-defined time intervals | TissueMiner |
| mqf_* functions | set of multi-query functions to quantify cell dynamics | TissueMiner |
The mqf functions constitute an important part of the TissueMiner API. They will be discribed thereafter.
openMovieDb provided with the TissueMiner API# Connection to the DB stored in the "db" variable
db <- openMovieDb(movieDir)
# Close the connection
dbDisconnect(db)
We can now overlay cells and vertices on the movie image. To do so, we built a dedicated render_frame() function that loads the specified frame of the time-lapse. This function takes the cell contour table and a desired frame as input variables. The render_frame() function alone returns the first layers of the graph that consists of a raster image of the wing and additional specifications such as the Y-axis flipping - scale_y_continuous(trans = “reverse”) - and the iso-scaling of the X and Y axes - coord_equal(). Additional options are appended using the + sign.
# Load the definitions of cell contours (calculated with the automated workflow `tm sm make_db`)
cellshapes <- locload(file.path(movieDir, "cellshapes.RData")) %>% print_head()
## Source: local data frame [6 x 5]
##
## frame cell_id x_pos y_pos bond_order
## (int) (int) (dbl) (dbl) (dbl)
## 1 0 10001 193 195 1
## 2 0 10001 189 199 2
## 3 0 10001 206 212 3
## 4 0 10001 229 206 4
## 5 0 10001 219 192 5
## 6 0 10001 217 190 6
## [1] 219714
# Plot cells and vertices on the original image
cellshapes %>%
# add overlay image (! connection to DB required !):
render_frame(50) +
# add a geometrical layer to represent cells as polygons
geom_polygon(aes(x_pos, y_pos, group=cell_id), color="green",fill=NA, size=0.2) +
# add a geometrical layer to represent vertices as points
geom_point(aes(x_pos, y_pos),color="red", size=0.4) +
# add a title to the graph
ggtitle("Cells and vertices overlaid on the image")
Other regions of interest can be manually defined by the user in Fiji (see Fiji macro). Of note, these additional ROIs are only taken into account if they were defined before running the automated workflow. Then force the automated workflow to rerun with the -R option, e.g tm sm -R roi_tracking that updates all previously performed calculations from the tracking of ROIs.
The automated workflow includes routines to browse the cell lineage and to follow ROIs in time once defined on a given image of the time-lapse. Please note that cells in contact with the margin are discarded because the segmentation and tracking quality isn’t optimum near the margin.
Example: Visualize cells in a selected ROI
# Load tracked ROIs:
lgRoiSmoothed <- locload(file.path(movieDir, "roi_bt/lgRoiSmoothed.RData")) %>%
print_head() %>%
filter(roi %in% c("cell_patch"))
## roi cell_id
## 1 cell_patch 10104
## 2 cell_patch 10110
## 3 cell_patch 10111
## 4 cell_patch 10112
## 5 cell_patch 10114
## 6 cell_patch 10116
## [1] 515
# Load cell shapes for plotting on the wing
cellshapes <- locload(file.path(movieDir, "cellshapes.RData"))
# Merge ROI with cell polygonal definition
cellshapesWithRoi <- cellshapes %>%
dt.merge(lgRoiSmoothed, by="cell_id", allow.cartesian=T) %>%
arrange(frame, cell_id, bond_order) ## .. because merge messes up the ordering
# Plot ROI on the wing
render_frame(cellshapesWithRoi, 50) +
geom_polygon(aes(x_pos, y_pos, fill=roi, group=cell_id), alpha=0.5) +
scale_fill_manual(values=c("cell_patch"="darkgreen"))
Videos are helpful to visualize the time evolution of patterns
Here, we use a parallelized loop over all frames of the time-lapse. The well-known avconv (formerly ffmpeg) program to create videos must be installed on your computer, see API requirements above.
# Make a video of the ROI on the wing
render_movie(cellshapesWithRoi, "bt_bhfix_peeled.mp4", list(
# geometrical layer to represent cells as polygons
geom_polygon(aes(x_pos, y_pos, fill=roi, group=cell_id), alpha=0.5),
# additional layer to fill cells with the green color
scale_fill_manual(values=c("cell_patch"="darkgreen"))))
TissueMiner provide a set of tools to quantify and visualize cell dynamics at different spatial scales. These tools are all prefixed with ‘mqf_’ as they perform multiple queries to the pre-processed data obtained with the TissueMiner automated workflow. Their common features are:
Mandatory argument: a path to a selected movie directory
Optional arguments: to control selected ROIs and other parameters that are specific to some subsets of mqf functions
These mqf_* functions are listed below. Just click on a function to get the detailed description of its options passed as arguments
| mqf_fg_* functions | Description |
|---|---|
| mqf_fg_nematics_cell_elong | get cell elongation nematics from the DB |
| mqf_fg_unit_nematics_CD | calculate cell division unit nematics |
| mqf_fg_unit_nematics_T1 | calculate cell neighbor change unit nematics |
| mqf_fg_cell_area | get cell area from the DB |
| mqf_fg_triangle_properties | get calculated triangle state properties |
| mqf_fg_bond_length | get bond length and positions from the DB |
| mqf_fg_cell_neighbor_count | calculate cell neighbor number from the DB |
| mqf_fg_dev_time | get developmental time from the configuration file |
| mqf_cg_roi_* functions | Description |
|---|---|
| mqf_cg_roi_cell_count | count cell number |
| mqf_cg_roi_cell_area | coarse-grain cell area |
| mqf_cg_roi_cell_neighbor_count | average cell neighbor count |
| mqf_cg_roi_polygon_class | average and trim cell polygon class |
| mqf_cg_roi_triangle_elong | coarse-grain cell elongation using triangles as a proxy |
| mqf_cg_roi_rate_CD | average cell division rate |
| mqf_cg_roi_rate_T2 | average extrusion rate |
| mqf_cg_roi_rate_T1 | average neighbor change rate |
| mqf_cg_roi_rate_isotropic_contrib | coarse-grain isotropic tissue deformation and its cellular contributions |
| mqf_cg_roi_rate_shear | coarse-grain anisotropic tissue deformation and its cellular contributions |
| mqf_cg_roi_nematics_cell_elong | coarse-grain cell elongation nematics by ROI |
| mqf_cg_roi_unit_nematics_CD | coarse-grain division unit nematics |
| mqf_cd_roi_unit_nematics_T1 | coarse-grain neighbor change unit nematics |
| mqf_cg_grid_* functions | Description |
|---|---|
| mqf_cg_grid_nematics_cell_elong | coarse-grain cell elongation nematics (DB) |
| mqf_cg_grid_unit_nematics_CD | coarse-grain division unit nematics |
| mqf_cg_grid_unit_nematics_T1 | coarse-grain neighbor change unit nematics |
To render cell bonds one must get bonds and their respective positions. Here, is an example in which different related tables must be joined together to pool the relevant data to be plotted:
For the sake of clarity, we built a dedicated mqf_fg_bond_length() function to get bond properties along with bond positions for plotting.
# we use the movieDir variable defined above
bond_with_vx <- mqf_fg_bond_length(movieDir, "raw") %>% print_head()
## movie frame cell_id bond_id vertex_id.2 vertex_id.1 x_pos.1 y_pos.1 x_pos.2 y_pos.2 bond_length roi time_sec timeInt_sec time_shift dev_time
## 1 demo 0 10000 53 1 51 28 343 18 339 11.65690 raw 0 291 54000 15
## 2 demo 0 10000 84 51 72 33 340 28 343 6.24264 raw 0 291 54000 15
## 3 demo 0 10000 128 72 100 45 344 33 340 13.65690 raw 0 291 54000 15
## 4 demo 0 10000 134 100 104 56 342 45 344 11.82840 raw 0 291 54000 15
## 5 demo 0 10000 135 104 105 60 337 56 342 6.65685 raw 0 291 54000 15
## 6 demo 0 10000 165 105 125 79 338 60 337 20.24260 raw 0 291 54000 15
## [1] 115219
bond_with_vx %>%
render_frame(50) +
# bonds are represented by segments using geom_segment
geom_segment(aes(x=x_pos.1, y=y_pos.1, xend=x_pos.2, yend=y_pos.2,
color=bond_length), # Here bond_length values are mapped to the color
size=0.5, lineend="round") +
# we overwrite the default color map by a custom rainbow palette
scale_color_gradientn(name="bond_length",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(bond_with_vx$bond_length, probs = 99/100)),
na.value = "red") +
ggtitle("Color-coded pattern of bond length")
# Here use the render_movie function
render_movie(bond_with_vx, "BondLengthPattern.mp4", list(
geom_segment(aes(x=x_pos.1, y=y_pos.1,xend=x_pos.2, yend=y_pos.2,color=bond_length),
size=0.3, lineend="round") ,
scale_color_gradientn(name="bond_length",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(bond_with_vx$bond_length, probs = 99/100)),
na.value = "red") # outliers are red
))
cellArea <- mqf_fg_cell_area(movieDir, rois = c("raw"), cellContour = T)
# Send a SQL query to get cell area in each frame
cellArea <- dbGetQuery(db,"select cell_id, frame, area from cells where cell_id!=10000") %>%
# add cell polygons into it
dt.merge(locload(file.path(movieDir, "cellshapes.RData")),
by = c("frame", "cell_id")) %>%
# order vertices around the cell contour for ploting cells as polygons
arrange(cell_id, frame, bond_order)
cellArea %>%
render_frame(50) +
# we now map cell area values to a color palette: fill=area
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area), alpha=0.7) +
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(cellArea$area, probs = 99.9/100)),
na.value = "red") +
ggtitle("Cell area pattern")
render_movie(cellArea, "CellAreaPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area)),
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
# ignore possible outliers by trimming the area distribution
limits=c(0,quantile(cellArea$area, probs = 99.9/100)),
na.value = "red")
))
# We now select the blade ROI that we defined using the draw_n_get_ROIcoord.ijm Fiji macro.
cellArea <- mqf_fg_cell_area(movieDir, rois = c("raw"), cellContour = T)
whole_tissue_roi <- locload(file.path(movieDir, "roi_bt/lgRoiSmoothed.RData")) %>%
filter(roi=="whole_tissue") %>% print_head()
## roi cell_id
## 1 whole_tissue 10000
## 2 whole_tissue 10102
## 3 whole_tissue 10104
## 4 whole_tissue 10107
## 5 whole_tissue 10109
## 6 whole_tissue 10110
## [1] 372
cellAreaInROI <- cellArea %>%
# A inner-join operation intersects the data
dt.merge(whole_tissue_roi, by = "cell_id")
cellAreaInROI %>%
render_frame(3) +
# we now map cell area values to a color palette: fill=area
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area)) +
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(cellAreaInROI$area, probs = 99.9/100)),
na.value = "red") +
ggtitle("Cell area pattern in whole_tissue")
cellShapesElong <- mqf_fg_nematics_cell_elong(movieDir, "raw", cellContour = T) %>% print_head()
## frame cell_id movie center_x center_y elong_xx elong_xy roi phi norm displayFactor x1 y1 x2 y2 time_sec
## 1 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## 2 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## 3 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## 4 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## 5 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## 6 0 10001 demo 209.731 201.607 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.1512 213.4106 202.0628 0
## timeInt_sec time_shift dev_time x_pos y_pos bond_order
## 1 291 54000 15 193 195 1
## 2 291 54000 15 189 199 2
## 3 291 54000 15 206 212 3
## 4 291 54000 15 229 206 4
## 5 291 54000 15 219 192 5
## 6 291 54000 15 217 190 6
## [1] 216795
cellShapesElong %>%
render_frame(50) +
# we now map cell elongation values to a color palette: fill=elongNorm
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=norm)) +
scale_fill_gradientn(name="elongation",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,1),
na.value = "red") +
ggtitle("Cell elongation pattern")
render_movie(cellShapesElong, "CellElongationPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=elongNorm)),
scale_fill_gradientn(name="elongation",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(cellShapesElong$elongNorm, probs = 99.9/100)),
na.value = "red")
))
Using the TissueMiner grammar, one can get cell elongation nematics for plotting in one step
cellElongNematics <- mqf_fg_nematics_cell_elong(movieDir, "raw") %>% print_head()
## movie frame cell_id center_x center_y elong_xx elong_xy roi phi norm displayFactor x1 y1 x2 y2 time_sec
## 1 demo 0 10001 209.731 201.6070 0.296725 0.0746626 raw 0.1232524 0.3059742 24.23546 206.0514 201.15117 213.4106 202.06283 0
## 2 demo 0 10002 577.155 247.0320 0.340883 -0.0941206 raw 6.1484872 0.3536381 24.23546 572.9085 247.60748 581.4015 246.45652 0
## 3 demo 0 10003 532.862 143.6060 0.248538 0.2278700 raw 0.3710212 0.3371882 24.23546 529.0541 142.12457 536.6699 145.08743 0
## 4 demo 0 10004 545.181 242.5890 0.391461 0.0907914 raw 0.1139501 0.4018517 24.23546 540.3430 242.03532 550.0190 243.14268 0
## 5 demo 0 10005 548.963 48.8402 0.498495 -0.1980080 raw 6.0941352 0.5363809 24.23546 542.5791 50.06167 555.3469 47.61873 0
## 6 demo 0 10006 160.281 206.3480 0.147981 0.0344489 raw 0.1143596 0.1519378 24.23546 158.4519 206.13791 162.1101 206.55809 0
## timeInt_sec time_shift dev_time
## 1 291 54000 15
## 2 291 54000 15
## 3 291 54000 15
## 4 291 54000 15
## 5 291 54000 15
## 6 291 54000 15
## [1] 36261
cellElongNematics %>%
# crop a the image by defining squareRoi
render_frame(50) +
# plot nematics as segments
geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),
size=1, alpha=0.7, lineend="round", color="red", na.rm=T) +
ggtitle("Cell elongation pattern")
TissueMiner provides a coarseGrid() function to map cell positions to each element of a square grid.
TissueMiner provides more specific routines using this coarseGrid() function to average nematics in space and time. Concerning cell elongtation, the mqf_cg_grid_nematics_cell_elong() TissueMiner function allows the user to directly get coarse-grained nematics that are automatically scaled for display on the original image. The display factor can also be manually defined.
cellElongNematicsCG <- mqf_cg_grid_nematics_cell_elong(movieDir, gridSize = 96) %>%
print_head()
## movie frame xGrid yGrid roi cgExx_smooth cgExy_smooth phi norm x1 y1 x2 y2 time_sec timeInt_sec time_shift
## 1 demo 0 145 145 raw 0.3348337 -0.03992370 6.22384823 0.3372054 102.1862 147.5434 187.8138 142.4566 0 291 54000
## 2 demo 0 145 241 raw 0.3274483 0.09651457 0.14331552 0.3413758 102.0255 234.7986 187.9745 247.2014 0 291 54000
## 3 demo 0 241 145 raw 0.3615719 -0.03191985 6.23915905 0.3629781 194.8775 147.0319 287.1225 142.9681 0 291 54000
## 4 demo 0 241 241 raw 0.3249297 0.04857413 0.07419615 0.3285403 199.3279 237.9024 282.6721 244.0976 0 291 54000
## 5 demo 0 337 145 raw 0.3248126 -0.02979818 6.23744354 0.3261766 295.5569 146.8970 378.4431 143.1030 0 291 54000
## 6 demo 0 337 241 raw 0.3633810 0.06902224 0.09385421 0.3698781 290.1622 236.5911 383.8378 245.4089 0 291 54000
## dev_time
## 1 15
## 2 15
## 3 15
## 4 15
## 5 15
## 6 15
## [1] 700
render_frame(cellElongNematicsCG, 1) +
geom_segment(aes(x=x1, y=y1,xend=x2, yend=y2),
size=2, lineend="round", color="red", na.rm=T)
render_movie(cellElongNematicsCG, "CellElongationNematicPattern.mp4", list(
geom_segment(aes(x=x1, y=y1, xend=x2, yend=y2),
size=2, lineend="round", color="red", na.rm=T)
))
Using the TissueMiner grammar, one can get the cell neighbor count for plotting in one step
cellPolygonClass <- mqf_fg_cell_neighbor_count(movieDir, "raw", cellContour = T)
# Define a discrete color palette of polygon class
polygonClassColors <- c("2"="black","3"="darkgrey", "4"="green",
"5"="yellow", "6"="grey", "7"="blue",
"8"="red", "9"="purple", ">9"="black")
cellPolygonClass %>%
render_frame(60) +
geom_polygon(aes(x_pos, y_pos, fill=as.character(polygon_class_trimmed),
group=cell_id), alpha=0.7) +
scale_fill_manual(name="polygon class", values=polygonClassColors, drop=F) +
ggtitle("PolygonClassPattern")
# Define a discrete color palette of polygon class
polygonClassColors <- c("2"="black","3"="darkgrey", "4"="green",
"5"="yellow", "6"="grey", "7"="blue",
"8"="red", "9"="purple", ">9"="black")
render_movie(cellPolygonClass, "CellPackingPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, fill=as.character(polygon_class_trimmed),
group=cell_id), alpha=0.7),
scale_fill_manual(name="polygon class", values=polygonClassColors, drop=F)
))
# limit generation to a more reasonable range
genColors =c("0"="white", "1" = "red", "2"="green", "3"="cyan", ">3"="magenta")
# Send a SQL query to get the cell lineage
cellsWithLin <- dbGetQuery(db, "select cell_id, lineage_group, generation
from cell_histories") %>%
# fix a generation cut off for a more reasonable range
mutate(generation_cutoff=ifelse(generation>3, ">3", ac(generation))) %>%
# add cell vertices for cell rendering as polygons
dt.merge(locload(file.path(movieDir, "cellshapes.RData")), by = "cell_id") %>%
arrange(frame, cell_id, bond_order) %>% print_head()
## cell_id lineage_group generation generation_cutoff frame x_pos y_pos bond_order
## 1 10001 lg_2 0 0 0 193 195 1
## 2 10001 lg_2 0 0 0 189 199 2
## 3 10001 lg_2 0 0 0 206 212 3
## 4 10001 lg_2 0 0 0 229 206 4
## 5 10001 lg_2 0 0 0 219 192 5
## 6 10001 lg_2 0 0 0 217 190 6
## [1] 219714
cellsWithLin %>%
render_frame(65) +
geom_polygon(aes(x_pos, y_pos, fill=as.factor(generation_cutoff), group=cell_id), alpha=0.5) +
scale_fill_manual(name="generation", values=genColors) +
ggtitle("CellDivisionGeneration")
# Define a discrete color palette of polygon class
render_movie(cellsWithLin, "CellGenerationPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, fill=as.factor(generation_cutoff),
group=cell_id), alpha=0.5),
scale_fill_manual(name="generation", values=genColors)
))
Using the TissueMiner grammar, one can get cell division unit nematics for plotting in one step
cdNematics <- mqf_fg_unit_nematics_CD(movieDir, rois = "raw", cellContour = T) %>%
print_head()
## frame cell_id movie mother_cell_id roi normCDxx normCDxy phi center_x center_y x1 y1 x2 y2 time_sec timeInt_sec
## 1 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## 2 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## 3 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## 4 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## 5 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## 6 1 10210 demo 10209 raw -0.136106 -0.9906943 5.429522 621.598 289.8545 605.6698 308.1206 637.5262 271.5884 291 290
## time_shift dev_time variable x_pos y_pos bond_order
## 1 54000 15.08083 left_daughter_cell_id 605 288 1
## 2 54000 15.08083 left_daughter_cell_id 602 296 2
## 3 54000 15.08083 left_daughter_cell_id 616 307 3
## 4 54000 15.08083 left_daughter_cell_id 624 304 4
## 5 54000 15.08083 left_daughter_cell_id 632 298 5
## 6 54000 15.08083 left_daughter_cell_id 629 292 6
## [1] 298
## Plot CD nematics on image #70
cdNematics %>%
render_frame(5) +
geom_polygon(aes(x_pos, y_pos, group=cell_id), alpha=0.5, fill="blue", color="grey") +
geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),
size=1, lineend="round", color="orange", na.rm=T) +
ggtitle("Cell division unit nematic")
Using the TissueMiner grammar, one can coarse-grain cell division unit nematics for plotting in one step
cgCDnematicsSmooth <- mqf_cg_grid_unit_nematics_CD(movieDir, gridSize = 96, kernSize = 5)
## Plot CD nematics on image #55
cgCDnematicsSmooth %>%
render_frame(5) +
geom_segment(aes(x=x1, y=y1, xend=x2, yend=y2),
size=2, alpha=0.7, lineend="round", color="orange", na.rm=T) +
ggtitle("Coarse-grained cell division nematics")
# Define a discrete color palette of polygon class
render_movie(cgCDnematicsSmooth, "cgCDnematics.mp4", list(
geom_segment(aes(x=x1, y=y1, xend=x2, yend=y2),
size=1, alpha=0.7, lineend="round", color="orange", na.rm=T)
))
# Load the detected changes in topology
topoChangeSummary <- locload(file.path(movieDir, "topochanges/topoChangeSummary.RData")) %>%
select(cell_id, frame, num_t1_gained, num_t1_lost, num_neighbors_t)
# Filter T1 transitions and bring cellshapes for plotting
csWithTopoT1 <- locload(file.path(movieDir, "cellshapes.RData")) %>%
dt.merge(topoChangeSummary, allow.cartesian=TRUE) %>%
filter(num_t1_gained>0 | num_t1_lost>0) %>%
mutate(t1_type=ifelse(num_t1_gained>0,
ifelse(num_t1_lost>0, "t1 gain and loss", "t1 gain"), "t1 loss")) %>%
print_head()
## frame cell_id x_pos y_pos bond_order num_t1_gained num_t1_lost num_neighbors_t t1_type
## 1 0 10014 500 263 1 0 1 7 t1 loss
## 2 0 10014 524 269 2 0 1 7 t1 loss
## 3 0 10014 530 263 3 0 1 7 t1 loss
## 4 0 10014 517 253 4 0 1 7 t1 loss
## 5 0 10014 510 251 5 0 1 7 t1 loss
## 6 0 10014 486 255 6 0 1 7 t1 loss
## [1] 28408
## define t1 color attribute
T1cols <- create_palette(unique(csWithTopoT1$t1_type))
csWithTopoT1 %>%
render_frame(40) +
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=t1_type), alpha=0.5) +
scale_fill_manual(values = T1cols, drop = FALSE) +
ggtitle("Cell neighbor exchanges")
# Define a discrete color palette of polygon class
render_movie(csWithTopoT1, "CellNeighborChangePattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=t1_type), alpha=0.5),
scale_fill_manual(values = T1cols, drop = FALSE)
))
Using the TissueMiner grammar, one can get T1 unit nematics for plotting in one step
t1nematics <- mqf_fg_unit_nematics_T1(movieDir, "raw", cellContour = T) %>% print_head()
## frame cell_id movie type roi unit_T1xx unit_T1xy phi center_x center_y x1 y1 x2 y2 time_sec timeInt_sec time_shift
## 1 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## 2 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## 3 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## 4 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## 5 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## 6 0 10014 demo loss raw 0.9998438 -0.01767349 6.274348 491.324 259.1675 479.2067 259.2746 503.4413 259.0604 0 291 54000
## dev_time variable t1_type x_pos y_pos bond_order
## 1 15 cell_id loss 500 263 1
## 2 15 cell_id loss 524 269 2
## 3 15 cell_id loss 530 263 3
## 4 15 cell_id loss 517 253 4
## 5 15 cell_id loss 510 251 5
## 6 15 cell_id loss 486 255 6
## [1] 30557
T1cols = c("gain"="red", "loss"="green", "gain_and_loss"="blue")
t1nematics %>%
render_frame(20) + #
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=t1_type), alpha=0.5, color="grey") +
scale_fill_manual(values = T1cols, drop = FALSE) +
geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),
size=1, alpha=0.7, lineend="round", color="red", na.rm=T) +
ggtitle("Cell neighbor exchanges")
Using the TissueMiner grammar, one can coarse-grain T1 unit nematics for plotting in one step
cgT1nematics <- mqf_cg_grid_unit_nematics_T1(movieDir, rois="raw", gridSize = 96, kernSize = 11) %>%
print_head()
## movie frame xGrid yGrid roi cgT1xx_smooth cgT1xy_smooth phi norm scaledFact x1 y1 x2 y2 time_sec timeInt_sec time_shift dev_time
## 1 demo 0 145 145 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## 2 demo 0 145 241 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## 3 demo 0 241 145 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## 4 demo 0 241 241 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## 5 demo 0 337 145 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## 6 demo 0 337 241 raw NA NA NA NA 121.8394 NA NA NA NA 0 291 54000 15
## [1] 695
## Plot CD nematics on image #50
cgT1nematics %>%
render_frame(50) +
geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),
size=2, alpha=0.7, lineend="round", color="red", na.rm=T) +
ggtitle("Coarse-grained T1 nematics")
# Define a discrete color palette of polygon class
render_movie(cgT1nematics, "cgT1nematics.mp4", list(
geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),
size=1, alpha=0.7, lineend="round", color="red", na.rm=T)
))
CAUTION: Tissue orientation matters for nematic components. Indeed, nematic tensors are symmetric traceless tensors that are characterized by 2 components projected onto the x and y axis of the Cartesian system. In order to compare nematics amongst different time-lapses one has to make sure that the tissues have a similar orientation with respect to the x and y axes. In the workflow, one has the possibility to rotate the images along with the data see Fiji macro to obtain visually comparable time-lapses.
CAUTION: Cumulative quantities are strongly influenced by the developmental time. Therefore, movies must be aligned in time prior to comparison between movies. We have aligned the three WT wing movies in time by aligning the peaks of their respective average cell elongation curves as a function of time. One movie is used as a reference and time shifts are applied to other movies. These time shifts must be stored in a configuration file containing the algnModel table as defined here.
# In your unix terminal (or Docker QuickStart Terminal), type in:
# Dataset WT_1 (~800Mb)
curl https://cloud.mpi-cbg.de/index.php/s/SgxxQk5CkIpTLPW/download | tar -zxvf -
# Dataset WT_2 (~800Mb)
curl https://cloud.mpi-cbg.de/index.php/s/Z6ZR1b0sGWnC8Cj/download | tar -zxvf -
# Dataset WT_3 (~600Mb)
curl https://cloud.mpi-cbg.de/index.php/s/4BJiyxnCS7HFyKB/download | tar -zxvf -
# These datasets are automatically extracted in the example_data folder where the demo was also downloaded.
# Define path to all processed movies: TO BE EDITED BY THE USER
movieDbBaseDir="/Users/retourna/example_data"
# Define a working directory where to save the analysis:
outDataBaseDir=file.path(movieDbBaseDir, "multi-movie_analysis")
scriptsDir=Sys.getenv("TM_HOME")
# Load TissueMiner libraries
source(file.path(scriptsDir, "commons/TMCommons.R"))
source(file.path(scriptsDir, "commons/BaseQueryFunctions.R"))
source(file.path(scriptsDir, "commons/TimeFunctions.R"))
# Force the use the flywing_tm_config.R otherwise another config file may be used by default
source(file.path(scriptsDir, "config/flywing_tm_config.R"))
# Set up working directory
mcdir(outDataBaseDir)
Principle:
Usage:
Arguments:
Returned value:
We can now use any previously described mqf_* function in combination with the multi_db_query() function to apply this mqf function to multiple movies. For example, one can get the individual cell area for the 3 datasets, in the ROI “raw” that corresponds to all tracked cells.
Next, we apply the synchronize_frames() function that will chunk the time in discrete time intervals and match the closest frame to be nearly centered on each time interval.
# Define the movies to be compared
movieDirs <- file.path(movieDbBaseDir, c("WT_1","WT_2","WT_3"))
# Get the area of each individual cell for the 3 selected datasets
syncMovies <- multi_db_query(movieDirs, mqf_fg_cell_area, rois="raw", cellContour=T) %>%
# chunk time into 900sec adjacent time windows,
# and for each movie, get the closest frame to be centered on each interval
synchronize_frames(900) %>% print_head()
## Source: local data frame [6 x 4]
## Groups: movie [1]
##
## movie interval_mid syncFrame dev_time
## (fctr) (dbl) (dbl) (dbl)
## 1 WT_1 53550 0 15.00000
## 2 WT_1 54450 2 15.15886
## 3 WT_1 55350 5 15.39155
## 4 WT_1 56250 8 15.62622
## 5 WT_1 57150 11 15.86092
## 6 WT_1 58050 14 16.09716
## [1] 200
## movie frame cell_id area center_x center_y roi time_sec timeInt_sec time_shift x_pos y_pos bond_order interval_mid dev_time max_interval
## 1 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 141.6083 527.6415 1 59850 16.63721 114750
## 2 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 148.8365 547.2311 2 59850 16.63721 114750
## 3 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 157.2551 553.7217 3 59850 16.63721 114750
## 4 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 185.3934 538.9490 4 59850 16.63721 114750
## 5 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 186.2675 536.8909 5 59850 16.63721 114750
## 6 WT_1 21 10001 867.5 162.3585 536.0844 raw 5893 278 54000 169.3064 521.9135 6 59850 16.63721 114750
## min_interval
## 1 53550
## 2 53550
## 3 53550
## 4 53550
## 5 53550
## 6 53550
## [1] 1602009
We now render the 3 movies corresponding to the 3 datasets.
# use zip option of render_movie() to combine movies in Fiji
movieDir <- file.path(movieDbBaseDir, c("WT_1"))
db <- openMovieDb(movieDir)
render_movie(syncMovies %>% filter(movie=="WT_1"),
"WT_1_CellAreaPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area)),
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(syncMovies$area, probs = 99.9/100)),
na.value = "red")
), createZip = T)
dbDisconnect(db)
movieDir <- file.path(movieDbBaseDir, c("WT_2"))
db <- openMovieDb(movieDir)
render_movie(syncMovies %>% filter(movie=="WT_2"),
"WT_2_CellAreaPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area)),
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(syncMovies$area, probs = 99.9/100)),
na.value = "red")
), createZip = T)
dbDisconnect(db)
movieDir <- file.path(movieDbBaseDir, c("WT_3"))
db <- openMovieDb(movieDir)
render_movie(syncMovies %>% filter(movie=="WT_3"),
"WT_3_CellAreaPattern.mp4", list(
geom_polygon(aes(x_pos, y_pos, group=cell_id, fill=area)),
scale_fill_gradientn(name="area (px)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(0,quantile(syncMovies$area, probs = 99.9/100)),
na.value = "red")
), createZip = T)
dbDisconnect(db)
For these 3 datasets, we manually defined 2 ROIs. Each type of ROI is comparable between wings:
# Define a list of movies to compare (paths to ùovie directories)
movieDirs <- file.path(movieDbBaseDir, c("WT_1",
"WT_2",
"WT_3"))
# Define ROIs to be analyzed
selectedRois=c("whole_tissue","interL2-L3", "distL3")
We now query multiple databases to extract the information about cell area averaged over the selected ROIs.
avgCellArea <- multi_db_query(movieDirs, mqf_cg_roi_cell_area, selectedRois) %>% print_head()
## movie frame roi area.avg area.sum nbcell time_sec timeInt_sec time_shift dev_time
## 1 WT_1 0 distL3 538.0098 27438.5 51 0 287 54000 15.00000
## 2 WT_1 0 whole_tissue 762.2676 227918.0 299 0 287 54000 15.00000
## 3 WT_1 1 distL3 552.2745 28166.0 51 287 286 54000 15.07972
## 4 WT_1 1 whole_tissue 756.8944 229339.0 303 287 286 54000 15.07972
## 5 WT_1 2 distL3 552.5098 28178.0 51 573 276 54000 15.15917
## 6 WT_1 2 whole_tissue 758.2180 231256.5 305 573 276 54000 15.15917
## [1] 1206
Next, we plot the average cell area as a function of time for each movie and each ROI on the same graph for comparison.
ggplot(avgCellArea, aes(dev_time, area.avg*(0.208^2), color=movie)) +
geom_line() +
ylab(expression(paste("<Cell area> [",mu,m^2,"]"))) +
scale_color_manual(values = movieColors) +
facet_wrap(~roi, ncol=4) +
ggtitle("averaged cell area")
We now query multiple databases to extract the information about cell elongation nematics averaged over the selected ROIs.
# Query the DBs and calculate the norm of the cell elongation nematics for each cell
avgCellElong <- multi_db_query(movieDirs, mqf_cg_roi_nematics_cell_elong, selectedRois) %>%
print_head()
## movie frame roi roi_center_x roi_center_y cgExx_smooth cgExy_smooth phi norm x1 y1 x2 y2 time_sec
## 1 WT_1 0 distL3 551.6417 902.6196 -0.09564322 0.05275300 1.318776 0.10922685 543.0212 869.1414 560.2622 936.0979 0
## 2 WT_1 0 whole_tissue 494.2077 885.0022 -0.05426890 0.05255105 1.186137 0.07554288 485.2359 862.8400 503.1795 907.1644 0
## 3 WT_1 1 distL3 554.4271 903.9248 -0.08984036 0.03116325 1.403853 0.09509174 549.4260 874.2467 559.4282 933.6029 287
## 4 WT_1 1 whole_tissue 496.2994 884.4222 -0.04893280 0.04743912 1.185846 0.06815342 488.1993 864.4303 504.3994 904.4142 287
## 5 WT_1 2 distL3 557.1707 905.6653 -0.08443210 0.03722514 1.363171 0.09227400 551.1505 877.0878 563.1909 934.2428 573
## 6 WT_1 2 whole_tissue 498.2760 885.0726 -0.04308360 0.04740089 1.154259 0.06405498 490.0734 866.5327 506.4785 903.6125 573
## timeInt_sec time_shift dev_time
## 1 287 54000 15.00000
## 2 287 54000 15.00000
## 3 286 54000 15.07972
## 4 286 54000 15.07972
## 5 276 54000 15.15917
## 6 276 54000 15.15917
## [1] 1206
Next, we plot the norm of the average cell elongation nematics as a function of time for each movie and each ROI on the same graph for comparison.
ggplot(avgCellElong, aes(dev_time, norm, color=movie)) +
geom_line() +
ylab(expression(paste("<cell elongation norm>"))) +
scale_color_manual(values = movieColors) +
facet_wrap(~roi, ncol=4) +
ggtitle("averaged cell elongation norm")
We now query multiple databases to extract the information about cell neighbor number averaged over the selected ROIs.
avgNeighborNb <- multi_db_query(movieDirs, mqf_cg_roi_cell_neighbor_count, selectedRois) %>%
print_head()
## movie frame roi avg_num_neighbors time_sec timeInt_sec time_shift dev_time
## 1 WT_1 0 distL3 5.941176 0 287 54000 15.00000
## 2 WT_1 0 whole_tissue 6.006689 0 287 54000 15.00000
## 3 WT_1 1 distL3 5.941176 287 286 54000 15.07972
## 4 WT_1 1 whole_tissue 6.006601 287 286 54000 15.07972
## 5 WT_1 2 distL3 5.941176 573 276 54000 15.15917
## 6 WT_1 2 whole_tissue 6.016393 573 276 54000 15.15917
## [1] 1200
Next, we plot the average cell neighbor number as a function of time for each movie and each ROI on the same graph for comparison.
ggplot(avgNeighborNb, aes(dev_time, avg_num_neighbors, color=movie)) +
geom_line() +
ylab(expression(paste("<cell neighbor number>"))) +
scale_color_manual(values = movieColors) +
facet_wrap(~roi, ncol=4) +
ggtitle("averaged cell neighbor number")
We now query multiple databases to extract the information about cell neighbor number for each individual cell for all selected ROIs.
avgPgClass <- multi_db_query(movieDirs, mqf_fg_cell_neighbor_count, selectedRois, polygon_class_limit=c(3,9)) %>% print_head()
## movie frame cell_id neighbor_number polygon_class_trimmed roi time_sec timeInt_sec time_shift dev_time
## 1 WT_1 0 10005 6 6 whole_tissue 0 287 54000 15
## 2 WT_1 0 10006 4 4 whole_tissue 0 287 54000 15
## 3 WT_1 0 10009 6 6 whole_tissue 0 287 54000 15
## 4 WT_1 0 10010 6 6 whole_tissue 0 287 54000 15
## 5 WT_1 0 10017 7 7 whole_tissue 0 287 54000 15
## 6 WT_1 0 10019 5 5 distL3 0 287 54000 15
## [1] 426461
Next, we plot the cell neighbor number distribution per ROI
ggplot(avgPgClass, aes(ac(polygon_class_trimmed), fill=as.factor(polygon_class_trimmed))) +
geom_bar(color="white") +
scale_fill_manual(values=c("3"="black", "4"="green", "5"="yellow", "6"="grey", "7"="blue", "8"="red","9"="purple"), name="polygon class") +
xlab("Cell neighbor number") +
facet_grid(movie~roi) +
ggtitle("cell_neighbor_number_distribution")
To plot the cell division rate, after querying the databases, we chunk the time into discrete intervals of 1hour and average the division rates within each interval to smooth the data.
CDrateByTimeIntervals <- multi_db_query(movieDirs, mqf_cg_roi_rate_CD, selectedRois) %>%
chunk_time_into_intervals(3600) %>%
group_by(movie, roi,interval_mid) %>%
summarise(avgCDrate=mean(cell_loss_rate),
semCD=se(cell_loss_rate),
time_sec=interval_mid[1],
dev_time=mean(dev_time))
Next, we plot the division rate normalized per cell
ggplot(CDrateByTimeIntervals, aes(dev_time, avgCDrate, color=movie)) +
geom_line()+
geom_point(size=1, color="black") +
geom_errorbar(aes(ymin=(avgCDrate-semCD), ymax=(avgCDrate+semCD)),
size=0.3, width=0.4, color="black") +
ylab(expression(paste("CD rate [", cell^-1, h^-1,"]"))) +
scale_color_manual(values = movieColors) +
facet_wrap(~roi) +
ggtitle("CD rate")
To plot the cell neighbor change rate, after querying the databases, we chunk the time into discrete intervals of 1hour and average these rates within each interval to smooth the data.
T1rate <- multi_db_query(movieDirs, mqf_cg_roi_rate_T1, selectedRois) %>%
chunk_time_into_intervals(deltaT = 3600) %>%
group_by(movie, roi,interval_mid) %>%
summarise(avgT1rate=mean(cell_topo_rate),
semT1=se(cell_topo_rate),
time_sec=interval_mid[1],
dev_time=mean(dev_time))
Next, we plot the neighbor change rate normalized per cell
ggplot(T1rate, aes(dev_time, avgT1rate, color=movie)) +
geom_line()+
geom_point(size=1, color="black") +
geom_errorbar(aes(ymin=(avgT1rate-semT1), ymax=(avgT1rate+semT1)),
size=0.3, width=0.4, color="black") +
facet_wrap(~roi) +
ylab(expression(paste("T1 rate [", cell^-1, h^-1,"]"))) +
ylim(c(0.2,2.3)) +
scale_color_manual(values = movieColors) +
facet_wrap(~roi) +
ggtitle("T1 rate")
We now query multiple databases to extract the information about cell neighbor change orientation averaged for the “whole_tissue” ROI.
selectedRois=c("whole_tissue")
# Get cell neighbor change nematics and align movie start
T1Nematics <- multi_db_query(movieDirs, mqf_cg_roi_unit_nematics_T1, selectedRois) %>%
align_movie_start(movieDirs) %>%
mutate(frame=frame-closestFrame) %>%
group_by(movie) %>%
mutate(maxnormByMovie=max(norm,na.rm=T)) %>%
group_by(movie,roi) %>%
mutate(maxnormByRoi=max(norm,na.rm=T)) %>% print_head()
## Source: local data frame [6 x 20]
## Groups: movie, roi [1]
##
## movie frame roi roi_center_x roi_center_y cgT1xx_smooth cgT1xy_smooth phi norm x1 y1 x2 y2 time_sec
## (fctr) (int) (chr) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl) (dbl) (int)
## 1 WT_1 0 whole_tissue 462.5906 853.0220 -0.025804106 -0.007694537 4.85728687 0.02692690 461.3601 861.4551 463.8212 844.5890 5615
## 2 WT_1 1 whole_tissue 519.5607 942.1266 0.022688506 0.003717711 0.08120772 0.02299108 512.3080 941.5363 526.8134 942.7168 5893
## 3 WT_1 2 whole_tissue 477.2932 827.5738 -0.008210321 0.024650594 0.94615397 0.02598194 472.4842 820.9033 482.1023 834.2443 6171
## 4 WT_1 3 whole_tissue 430.5783 817.7019 0.017046864 0.006263175 0.17604958 0.01816103 424.9192 816.6952 436.2374 818.7086 6448
## 5 WT_1 4 whole_tissue 495.1280 945.1814 -0.050444578 0.025106067 1.33989495 0.05634687 491.0466 927.8209 499.2093 962.5419 6729
## 6 WT_1 5 whole_tissue 466.7140 786.9665 -0.046644580 0.001571496 1.55395726 0.04667104 466.4653 772.1973 466.9627 801.7358 7022
## Variables not shown: timeInt_sec (int), time_shift (dbl), dev_time (dbl), closestFrame (int), maxnormByMovie (dbl), maxnormByRoi (dbl)
## [1] 578
T1Nematics$roi <- factor(T1Nematics$roi, levels=c("whole_tissue","distL3"))
We render the average cell neighbor change orientation in a circular diagram. A similar plot can be done for the average cell division orientation (see the QuickStart tutorials).
ggplot(T1Nematics , aes()) +
geom_segment(aes(x=phi, y=0, xend=phi, yend=(norm), color=dev_time),size=1, alpha=0.5) +
geom_segment(aes(x=mod2pi(phi+pi), y=0, xend=mod2pi(phi+pi), yend=(norm),
color=dev_time), size=1, alpha=0.5) +
scale_color_gradientn(name="time (hAPF)",
colours=c("black", "blue", "green", "yellow", "red"),
limits=c(min(T1Nematics$dev_time),max(T1Nematics$dev_time)),
na.value = "red") +
coord_polar(start=-pi/2,direction=+1)+
scale_x_continuous(breaks=seq(0,3*pi/2,pi/2),
labels=c(expression(pi),expression(paste(pi/2," Ant")),
expression(0),expression(-pi/2)),
limits=c(0,2*pi)) +
xlab("") +
ylab("T1 nematic norm") +
facet_grid(movie~roi) +
theme(
plot.title = element_blank(),
panel.grid.minor=element_blank(),
plot.margin = unit(c(0,0,0,0), "cm")) +
ggtitle("T1 nematics - avg by roi in frame")
To optimally compare the isotropic deformation curves between movies, one can interpolate the data in time instead of relying on the supposedly constant frame rate of the microscope.
movieDirs <- file.path(movieDbBaseDir, c("WT_1","WT_2", "WT_3"))
selectedRois=c("whole_tissue", "distL3", "distInterL3-L4")
# query multiple databases
isoContrib <- multi_db_query(movieDirs, mqf_cg_roi_rate_isotropic_contrib, selectedRois) %>%
filter(isoContrib!="tissue_area")
# interpolate data
deltaT=60 # sampling (60 seconds)
avgIsoDefRateInterpolated <- isoContrib %>%
align_movie_start(movieDirs) %>%
mutate(min_dev_time=min(dev_time),
max_dev_time=max(dev_time)) %>%
group_by(movie, roi, isoContrib) %>%
# here, we do a linear interpolation of the data (60 sec per point)
do({
with(., approx(dev_time, value.ma,
xout = seq(min_dev_time[1], max_dev_time[1],
by = deltaT/3600), method = "linear")) %>% as.df()
}) %>%
rename(dev_time=x, value.ma=y) %>%
filter(!is.na(value.ma)) %>%
ungroup() %>% mutate(movieNb=length(unique(movie))) %>%
group_by(roi, isoContrib, dev_time) %>%
# Make sure that further calculations will be done on values present in all compared movies
filter(n()==movieNb) %>% ungroup()
# average data between the 3 movies amongst ROI
avgIsoDefRateSummary <- avgIsoDefRateInterpolated %>%
group_by(roi, isoContrib, dev_time) %>%
summarise(value.avg=mean(value.ma), value.sd=sd(value.ma))
We now plot the isotropic deformation rate averaged in each ROI and further averaged between movies. We also plot the standard deviation between movies.
# Plot average of iso contribution rates in 3 WT and their respective standard deviation
ggplot(avgIsoDefRateSummary, aes(dev_time, value.avg, color=isoContrib)) +
geom_line() +
geom_ribbon(aes(ymin=(value.avg-value.sd), ymax=(value.avg+value.sd), fill=isoContrib),
alpha=0.2, linetype="dotted", size=0.2) +
xlab("Time [hAPF]")+
scale_x_continuous(breaks=seq(16,36, 4),limits=c(16,32)) +
ylab(expression(paste("rate [",h^-1,"]"))) +
scale_color_manual(values=isotropColors) +
scale_fill_manual(values=isotropColors) +
facet_wrap(~roi) +
ggtitle("averaged rates of isotropic deformation")
We next calculate the cumulative isotropic deformation that we further average between movies.
avgIsoDefCum <- avgIsoDefRateInterpolated %>%
group_by(movie, roi, isoContrib) %>%
mutate(timeInt=c(0,diff(dev_time)), value.cumsum=cumsum(value.ma*timeInt)) %>%
group_by(roi, isoContrib, dev_time) %>%
summarise(cumsum.avg=mean(value.cumsum), cumsum.sd=sd(value.cumsum))
We now plot the cumulative isotropic deformation further averaged between movies. We also plot the standard deviation between movies.
CAUTION: movies must be well registered in time to get an optimal comparison between movies
ggplot(avgIsoDefCum, aes(dev_time, cumsum.avg, color=isoContrib)) +
geom_line() +
geom_ribbon(aes(ymin=(cumsum.avg-cumsum.sd), ymax=(cumsum.avg+cumsum.sd), fill=isoContrib),
alpha=0.2, linetype="dotted", size=0.2) +
xlab("Time [hAPF]")+
ylab("cumulative sum") +
scale_x_continuous(breaks=seq(16,36, 4),limits=c(16,32)) +
scale_y_continuous(breaks=seq(-2, 1, 0.4), limit=c(-1.2, 1)) +
scale_color_manual(values=isotropColors) +
scale_fill_manual(values=isotropColors) +
facet_wrap(~roi) +
ggtitle("cumulative avg isotropic deformation")
The changes in aspect ratio of cells and tissue correspond the pure shear that we term shear. The tissue decomposition into each type of cellular contributions is obtained by the triangulation method that we published elsewhere (Etournay et al., 2015). Here, we show the triangulation of the cell network.
movieDir <- file.path(movieDbBaseDir, c("WT_1"))
triProperties <- mqf_fg_triangle_properties(movieDir, "raw", triContour = T) %>%
print_head()
## tri_id movie frame ta_xx ta_xy ta_yx ta_yy tri_area s_a theta_a two_phi_a Q_a roi time_sec timeInt_sec time_shift
## 1 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 2 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 3 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 4 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## 5 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## 6 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## dev_time cell_id x_pos y_pos
## 1 18.22694 10001 104.59096 514.2104
## 2 18.22694 10008 73.94130 534.5956
## 3 18.22694 10179 67.50174 515.1043
## 4 15.00000 10001 176.60837 576.1622
## 5 15.00000 10008 158.56598 612.6568
## 6 15.00000 10320 152.27188 586.8624
## [1] 1680423
## tri_id movie frame ta_xx ta_xy ta_yx ta_yy tri_area s_a theta_a two_phi_a Q_a roi time_sec timeInt_sec time_shift
## 1 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 2 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 3 1 WT_1 41 25.73521 4.237465 -8.08430 12.82599 364.3366 2.949039 -0.3092837 5.684289 0.3459120 raw 11617 283 54000
## 4 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## 5 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## 6 2 WT_1 0 16.10049 4.141756 -17.93016 16.97369 347.5472 2.925450 -0.5884726 4.060672 0.3625545 raw 0 287 54000
## dev_time cell_id x_pos y_pos
## 1 18.22694 10001 104.59096 514.2104
## 2 18.22694 10008 73.94130 534.5956
## 3 18.22694 10179 67.50174 515.1043
## 4 15.00000 10001 176.60837 576.1622
## 5 15.00000 10008 158.56598 612.6568
## 6 15.00000 10320 152.27188 586.8624
## [1] 1680423
db <- openMovieDb(movieDir)
triProperties %>%
render_frame(40) +
geom_polygon(aes(x_pos, y_pos, group=tri_id),
fill="yellow", color="black", alpha=0.9, size=0.1)
To optimally compare the pure shear deformation curves between movies, one can interpolate the data in time instead of relying on the supposedly constant frame rate of the microscope.
# query multiple databases
shearData <- multi_db_query(movieDirs, mqf_cg_roi_rate_shear, selectedRois)
shearRateSlim <- subset(shearData, (tensor=="CEwithCT" | tensor=="correlationEffects" |
tensor=="nu" | tensor=="ShearT1" |
tensor=="ShearT2" | tensor=="ShearCD"))
shearRateSlim$tensor <- factor(shearRateSlim$tensor,
levels=c("ShearCD", "CEwithCT", "correlationEffects",
"nu", "ShearT1", "ShearT2"),
labels=c("cell_division", "cell_elongation_change",
"correlation_effects","total_shear","T1", "T2"))
# interpolate data
deltaT=60
shearRateInterpolated <- shearRateSlim %>%
align_movie_start(movieDirs) %>%
select(-c(xy,yx,yy,xy.ma,yx.ma,yy.ma, TimeInt.ma,
phi, norm,time_sec,timeInt_sec,closestFrame,time_shift)) %>%
mutate(min_dev_time=min(dev_time),
max_dev_time=max(dev_time)) %>%
group_by(movie, roi, tensor) %>%
do({
with(., approx(dev_time, xx.ma,
xout = seq(min_dev_time[1], max_dev_time[1],
by = deltaT/3600), method = "linear")) %>% as.df()
}) %>%
rename(dev_time=x, XX=y) %>%
filter(!is.na(XX)) %>%
ungroup() %>% mutate(movieNb=length(unique(movie))) %>%
group_by(roi, tensor, dev_time) %>%
# Make sure that further calculations will be done on values present in all compared movies
filter(n()==movieNb) %>% ungroup()
# average data between the 3 movies amongst ROI
shearRateSummary <- shearRateInterpolated %>%
group_by(roi, tensor, dev_time) %>%
summarise(xx.avg=mean(XX), xx.sd=sd(XX))
We now plot the pure shear deformation averaged in each ROI and further averaged between movies. We also plot the standard deviation between movies.
# Plot avg and standard deviation for each tensor among 3 WT
ggplot(shearRateSummary, aes(dev_time,xx.avg*100, color=tensor)) +
geom_line() +
geom_ribbon(aes(ymin=100*(xx.avg-xx.sd), ymax=100*(xx.avg+xx.sd), fill=tensor),
alpha=0.2, linetype="dotted", size=0.2) +
xlab("Time [hAPF]")+
scale_x_continuous(breaks=seq(16,36, 2),limits=c(16,34)) +
ylab(expression(paste("shear rate xx [",10^-2,h^-1,"]"))) +
scale_color_manual(values=shearColors) +
scale_fill_manual(values=shearColors) +
facet_wrap(~roi) +
ggtitle("shear decomposition")
We next calculate the cumulative pure shear deformation that we further average between movies.
shearCumSumSummary <- shearRateInterpolated %>%
group_by(movie, roi, tensor) %>%
mutate(timeInt=c(0,diff(dev_time)),cumSum_xx=cumsum(XX*timeInt)) %>%
group_by(roi, tensor, dev_time) %>%
summarise(xxCumSum.avg=mean(cumSum_xx, na.rm = F), xxCumSum.sd=sd(cumSum_xx, na.rm = F))
We now plot the cumulative pure shear deformation further averaged between movies. We also plot the standard deviation between movies.
CAUTION: movies must be well registered in time to get an optimal comparison between movies
ggplot(shearCumSumSummary, aes(dev_time,xxCumSum.avg, color=tensor)) +
geom_line() +
geom_ribbon(aes(ymin=(xxCumSum.avg-xxCumSum.sd), ymax=(xxCumSum.avg+xxCumSum.sd), fill=tensor),
alpha=0.2, linetype="dotted", size=0.2) +
xlab("Time [hAPF]")+
scale_x_continuous(breaks=seq(16,36, 2),limits=c(16,34)) +
ylab(expression(paste("cumulative PD shear"))) +
scale_color_manual(values=shearColors) +
scale_fill_manual(values=shearColors) +
facet_wrap(~roi) +
ggtitle("cumulative shear decomposition")